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Motivated by various problems in physics and applied mathematics, we look for constraints and 
properties of real Fourier-positive functions, i.e. with positive Fourier transforms. Properties of 
the “Dirac comb” distribution and of its tensor products in higher dimensions lead to Poisson 
resummation, allowing for a useful approximation formula of a Fourier transform in terms of a limited 
number of terms. A connection with the Bochner theorem on positive definiteness of Fourier-positive 
functions is discussed. As a practical application, we find simple and rapid analytic algorithms for 
checking Fourier-positivity in 1- and (radial) 2-dimensions among a large variety of real positive 
functions. This may provide a step towards a classification of positive positive-definite functions. 

PACS numbers: 


I. INTRODUCTION 

We call “Fourier-positivity” a property of a real positive function whose Fourier transform is itself positive. Besides 
of a purely mathematical interest [1], the study of such pairs of functions is motivated by applications in physics 
and applied mathematics [2]. In physics, typically, both functions correspond to observables, i.e. measurable a priori 
positive quantities. Well-known examples exist in one- and especially two-dimensional cases. Let us quote for instance 
the Fourier-Bessel transform, i.e. radial version of the 2-dimensional Fourier transform, relating the gluon and dipole 
distributions [3] inside a hadron in the framework of Quantum Chromodynamics of strong particle interactions. 

There exists a fundamental property characterizing Fourier positivity, which uses the Bochner theorem [4]: “Fourier- 
positivity” of a real function tplf) is equivalent to the statement that ip is not only positive but also positive-definite. 
“Positive-definiteness” means that for any set of positions, {r^, f = 1,..., n}, the nxn matrix M with elements 
is positive definite, i.e. 


Ui ip{fi — fj) Uj > 0 Vu, Vn G N . (1) 

i j = i 

In other terms, the lowest eigenvalue of the matrix M remains positive for all u and all values of n. 

The Bochner theorem with its applications appears to be still the major tool in the domain. However, to our 
knowledge, there does not yet exist a mathematical classification of Fourier-positive functions which, for instance, 
could allow for an appropriate parametrization for model building. Testing positive-definiteness (I) cannot be done 
concretely, due to the generality of the constraints. Gonversely, numerically computing Fourier transforms for checking 
Fourier-positivity is obviously possible, but it does not give general or analytical, easy means to select a priori 
appropriate Fourier-positive sets of functions. 

Our approach is to find new constraints of Fourier-positivity allowing for simple and efficient selection rules of 
functions with positive Fourier transform ip, given the positive input ip. We thus consider a pair of real even functions 
on d-dimensional real vector spaces, tp{f) and which are taken to be d-dimensional Fourier transforms one from 
the other. 
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We have already performed preliminary studies on this problem. In our initial paper [5], we examined the distri¬ 
bution of Fourier-positive functions among arbitrary combinations of a finite basis of Fourier eigenfunctions using the 
algebra of Hermite polynomials (in one dimension) and the algebra of Laguerre polynomials (in radial two dimensions). 
The main outcome of this first study, using the Sturm algorithm on the number of polynomial zeros, is to reveal the 
rather intricate geometry of the manifold of solutions. In a second study [6], we derived generalized sufficient proper¬ 
ties, based on an extension of convexity conditions by analytic continuation of ip into the complex plane and Jensen 
inequalities. However, the set of obtained constraints proved to be too weak to reliably check Fourier-positivity in 
the testing domain. 

In the present work we propose a method providing a satisfactory detection of Fourier positivity, thoroughly tested 
for large sets of functions in one- and radial two dimensions. It is based on a different approach from aforementioned 
studies, using remarkable properties of “Dirac comb” mathematical distributions through Fourier transforms in any 
dimension. 

The plan of this paper is the following. In section II, we recall the definition and properties of the “Dirac comb” 
distribution and of its tensor products, leading to the well-known Poisson resummation formula. In section III, we 
then derive a useful approximation formula of a Fourier transform in term of a finite (and limited) sum of data on the 
candidate Fourier-positive function. This is shown, in section IV, to be equivalent to the application of the Bochner 
theorem in a finite regular lattice of points of varying size. In section V we apply these results to large random sets 
of functions, either Fourier-positive or not and show the relevance of our method to select Fourier-positivity with an 
analytic evaluation of the Poisson summation. The final section VI is devoted to a summary of results and prospects 
for a deeper understanding of Fourier-positivity. 


II. DIRAC COMBS AND THE POISSON RESUMMATION FORMULA 


As we shall see, a key ingredient of our approach is to take advantage of the Poisson resummation formulas, which 
can be easily derived from the so-called “Dirac comb” mathematical distribution. 


Y.d{r-k) = 

feez fcez 


fceZ'* 


fceZ'* 


( 3 ) 


where the second line is its d-dimensional tensor product. 

“Dirac combs” are formally invariant under Fourier transform, and also Fourier-positive. Indeed, as easily derived 
from the definition (2) inserted into the second line of (3) one writes. 


ydre— ^ = (27r)-^ E = E ' (4) 

feez<i kGZ‘‘ fceZ'^ 

The connection between Dirac combs and Fourier-positive functions is obtained by considering the “characteristic 
function” defined by: 

F{0,r) = Y ip{kiri,...,kdrd) . (5) 

feeZ'i 

As we shall discuss in further sections the condition for ip to be Fourier-positive transfers the condition to a positivity 
of F{9, r), namely. 


F(9,r) >0 Vr,0 G (g) [0, 27r['^ . 


In fact, noting from the second line of (2) that, 


V’(fciri,..., kdVd) = 


1 

(27r)‘^/^ 


ds e 




<f(s) , 


( 6 ) 


( 7 ) 
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we are able to use the rf-dimensional “Dirac comb” relation (3) in order to rewrite the positivity condition (6) as, 


F(0,f) = 


1 


(27r)‘^/2 

1 

(27r)‘^/^ 

(27r)‘^/2 

n.lilD-l 


ds if{s) 




feez‘d 


dsipis) \ Y 

, feGZ'* 


D . 

fC-j 


27r 


EH 


/ 27rfci + 6*1 27rA:d + Od 


ri 


Td 


> 0 Vf,d G [0,27r[^ . 


( 8 ) 


The equality of the two expressions (5) and (8) of F{d, r) is nothing but a version of the d-dimensional Poisson 
resummation formula [7]. 

An interesting insight on the properties of (8) is obtained by a change of variables {Oj,rj) {sj,rj) with 


0 -, 


Si = 


j = l,...,d; F{0,r) F{s,r) 


(9) 


The positivity condition (6) may thus be rewritten and renormalized in such a way as to read, 

F(.>1 s + .“ + E.>0. (10) 




fceZ'i 


Under this form, the Poisson resummation formula (10) allows for an interesting Fourier transform relation which can 
be qualitatively (and made quantitative in the next section) outlined as follows, 

• Approximation of the right-hand side of Eq.(lO). Assume that the function tpiri,... ,rd) has a hnite range R 
for each of its arguments, namely that it is negligible^ if any \rj\ > R. Then one can choose a positive integer K 
so that the right-hand summation can be truncated into a finite number of terms. Indeed, define a parameter, 
^’min = R/K, and consider only situations where \rj \ > rniinj Vj = 1,..., d. Clearly, every \kj\ becomes bounded 
by 77, hence. 


F{s,r) 


ID 


■ rd\ 


(27r)^/2 


Y ifikiri,... ,kdrd) 

feeZ'* ; \kj\<K,Vj 


( 11 ) 


Given K, there appears a minimal value of each rj for practical calculations if K remains fixed. 

• Approximation of the left-hand side of Eq.( 10). Assume also that the function tp{s) has a finite range for each of 
its arguments, |sj| < S', Vj, beyond which it is negligible. Then one finds that the left-hand summation can be 
limited to its first term /ij = 0 , j = 1,..., d. This only remaining significant term is just the Fourier transform 
of if, hence. 


E ‘dfsiT 


27rhi 


ri 


) SdF 


inhd 

Td 


— ‘P {.Sli ■ ■ ■ 1 Sd) ) 


( 12 ) 


This holds for some maximal value of each \rj\, depending on S. This maximum, rmax, will be derived in the 
next section. 

Let us comment these two approximations. Eq.(ll) is meant to obtain a good approximation of the characteristic 
function itself. Obviously, this expression exhibits an approximation formula for a Fourier transform, and, usually, a 
convergent result at the limit, K —>■ oo, or, as well, rmin —t 0.. This requires a priori K to be large enough. However, 
our aim is to look for a good enough approximation for a limited number of terms in (11). In this case the variable 
f'min plays the role of a “resolution” parameter, Ar = R/K, on the function if allowing to test its Fourier-positivity. 


^ It is assumed to be small enough even in a summation like in (10). 
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The approximation described by Eq. (12) is of different nature and is directly related to the properties of the Poisson 
resummation. Indeed, the trade of variables kiVi into Si + , typical of the Poisson resummation formula (10), 

allows for a rapid decrease rate in the left-hand series, for small enough values of 

By combining both approximations (11,12), one has the interesting approximation property, valid in a restricted 
domain for r, of a Fourier transform by a finite sum. 


</5(si,.--,S<i) 


Vi ■■■rd\ 


K 

kj—1 


kjTjSj 


(13) 


allowing to check the positivity and other properties of (p from a finite number of values of ip for the variables rj in a 
given range, rmin < |?"il < ?’max, delimited by both a lower bound (related to the right-hand side approximation) and 
an upper bound (related to the left-hand side approximation) induced from the Poisson resummation (10). 


III. FOURIER TRANSFORM via POISSON RESUMMATION 
A. The one-dimensional Fourier case 


Here we consider a conjugate pair of real even functions ipiYiVY) of variables which are Fourier-conjugated 
one with the other, see formulas (2) when d=l. We look for the application of the characteristic function F{s, r) and 
the corresponding Poisson resummation formula (10) to the one-dimensional problem of Fourier positivity, using the 
approximation scheme (13) outlined in the previous section. 

Let us first write the Poisson resummation formula (10) in the one-dimensional case. 


F{s,r) 


= E 




2Trh 

r 



Y, Yikr) . 

k^Tj 


(14) 


Under this form, the approximation properties leading to (13) can be made quantitative as follows; 


• Approximation of the right-hand side summation in EqflA). Let us consider a fixed value, moderately large, 
of the positive integer parameter K. Consider a function tpY) with a finite and non-zero range R, define the 
parameter, rmin = R/K^ and restrict r to be larger than rmin- This gives, 


F{s,r) 



E 

\k\<K 


^ikrs 


(15) 


This summation in (14) is thus limited^ to just terms with |fc| < K and looks like a discretized approximation, 
with step r, of the Fourier integral. 

• Approximation of the left-hand side summation in Eq.{14). Let us assume a finite range S of the function (/3(s), 
Fourier partner of namely (p{s) is negligible if |s| > S. 

Besides an obvious condition, |s| < S, the condition for retaining the first term, h = 0, as the only significant 
term in the left-hand side summation (14), is. 


s ± 


27r 




(16) 


since all other terms, with \h\ > 1, can then be neglected (provided one has good convergence properties). A 
bound on |r| is obtained as follows. Clearly, whatever the signs of r and s, the condition (16) reduces to, 


27r 




(17) 


^ We assume that the cut-off R is strong enough to ensure a fast convergence of the series (15). 
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with two branches, depending on the relative values of |s| and 27r/|r|, 

Ott ^tt 

-N + ^> 5 . ( 18 ) 

Since |s| < S, the hrst branch is useless. The second one gives 27r/|r| > 2S, because |s| may reach S. Accordingly, 
one finds the bound, |r| < rmax = tt/S'. 

By combining both approximations, one obtains the approximation property of a Fourier transform by a finite 
sum, namely, 


K 


ikrs 


‘p{s) ^ —f= e 

For (19) to be valid, one has to choose an appropriate range for |r|, namely, 

R I I tt 
K ^ ^ S ' 

This, in turn, requires the following condition on the truncation parameter, AT, 

RS 


K > 


TT 


(19) 


( 20 ) 


( 21 ) 


the value of which is not too large for a pair of Fourier partners whose cut-offs are such that the product, RS, 
is a finite number. 


B. The radial two-dimensional Fourier case 


The radial two-dimensional problem considers 'tjj{r) and a pair of radial conjugate functions on R’*', namely 


(•-1-00 


(p{s) = t;- I dr 'ip{r) = / rdr jQ(sr) ij){r). 


27r 

1 


(•-1-00 


7 /’( 7 ’) = TT I ds e (p{s) = / sds Jo{rs) if{s) 

2 .TT . 02 ./n 


( 22 ) 

(23) 


where r = and s = |s| denote here the radial variables of vectors in both conjugated 2-dimensional spaces. 

Following the general formalism of section 11, the two-dimensional Poisson resummation ensures the positivity of 
the 2-dimensional characteristic function for a positive Fourier transform namely, 


2 _ 

F{si,S2,r) = ^ (t-v/-I- ^t{msi+nS2)r ^ Q VrG[ 0 ,Oo[, Vsi,S2 • 


(24) 


m,nGZ 


As we discuss now, the condition (24) happens to furnish also a condition for Fourier-positivity of the function ipi''’) 
in the range where A(si,S 2 ;r) gives a direct link to its Fourier transform similar to the one-dimensional case. This 
comes again from a 2-dimensional Poisson summation formula relating 'tjj{r) to (p{s), namely. 


A(si,S 2 ,r) = Y \\l + S 2 )^ I > 0. 

hiM&z yv r r y 


(25) 


• Approximation of the summation in Eq.{2A). Let us consider a fixed value of the variable r. Then, we consider 
functions if with a finite range R. It is easy to see that the summation (24) can be limited^ to just terms with 
\m\, \n\ < K = [R/r] -|- 1 where [R/r] is the integer part of i?/r . 


2 _ 

F{si,S2,r) ~ Y V' (ry'm^ + Kr > R . 


(26) 


m,n^Z; |m|,|n|<X 


^ We assume that the cut-off R is strong enough to ensure a fast convergence of the series(24). 
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• Approximation of the summation in Eq.(25). Consider the case of a finite range S of the function (p, Fourier 
transform of ip. Namely, assume that ip{s) is negligible for s = \/s\ + s\ > S. Then, consider the set of points, 
{si, 52 }, inside the circle with radius S', namely, sf + < S. The conditions for retaining only the term with, 

hi = h 2 = 0, in (25) read, 



2 

+ ( 52 )^ > S^ (for hi = ±1, /i 2 


0) and (si)^ + 



> S^ (for h 2 = ±1, hi 


0) . (27) 


This is obtained if, 27r/r > 2S. Indeed, under such a condition for 27r/r, the points, {±27r/r + si, S 2 }, are pushed 
out of the circle which confines {si, 52 }- And the same holds for the points {si, ±27r/r + 52 }- Accordingly, the 
summation (25) can be reduced^ to its simplest term, h = 0, i.e., the only remaining significant term. It is just 
the Fourier transform ip{s) of ip. 

By combining both approximations, one obtains the approximation property of a Fourier transform by a finite 
sum, namely, 

^ (r\/m2 + n2) ^^(ms^+ns,)r (28) 

m,n€Z; |m|,|n|<K’ 

For (28) to be valid, one has to choose appropriate bounds for r, which are similar to those of the Id case, 
namely. 


R 

K 


< r < 


S 


K > 


RS 

TT 


(29) 


IV. FROM POISSON FORMULA TO BOCHNER POSITIVE-DEFINITENESS 

Our aim in this section is to exhibit the connection of the positivity condition (6) of the characteristic function 
F{0,r) with the positive-definiteness of certain moment matrices [8] related to the Bochner theorem [4]. This in 
turn implies the positivity of the lowest eigenvalue of these matrices (and thus also of the corresponding matrix 
determinants). For practical reasons, we shall focus the discussion on the one- and radial two-dimensional cases, but 
the method is general. 


The one-dimensional case 


Let us recall the one-dimensional Poisson formula (14) and its positivity condition under the form. 


A(6»,r) = '^ippkr) exp (i/cd) = ^ 




feez 


fcez 


2Trk + 9 


> 0, V (r,d) G (] — oo,oo[(g) [0,27r[) . 


(30) 


As discussed in the preceding section, the positivity condition (30) happens to be an equivalent formulation of the 
Fourier-positivity of the input function ip{r). In the following we show that this positivity condition can be rephrased 
in terms of a specific application of the Bochner theorem [4]. 

In the present case, starting from Eq.(30), the problem^ may be formulated as finding the conditions on the set of 
functions. 


iPikr) = 



F{9,r) dO , 


fc = 0,±l,±2,... 


(32) 


We assume again that the cut-off S is strong enough to ensure a fast convergence. 

^ In the mathematics literature, one may refer to the problem of moments [8, 9], and more specifically in our case to the trigonometric 
moment problem [8, 10, 11]. The trigonometric moment problem can be expressed as follows: Find a bounded, positive function 
F{0) > 0, 0 E [0, 27r] such that its trigonometric moments, 

Pn = — / e*"® F{6) dO , n = 0, ±1, ±2,..., p-n = TU Vn . (31) 

Stt Jq 


have a prescribed set of values. 



7 


such that the function F{0,r) be positive. 

Let us recall® a simple derivation of these conditions. Consider a real positive polynomial P{z,z), with 2 : being a 
complex variable, this polynomial being obtained as the squared modulus of an arbitrary complex polynomial Q{z), 

k 

P{z,z) = \Q{z)\‘^ = z^z’’ Cl > 0, V complex vector c . (34) 

j,i 

Then, choosing 2 ; = in Eq.(34) and integrating over F{6,r) as in Eq.(32), one finds, 

-I ^277 ^ ^ 

= — de P{F\e-^% Cl F{e,r) = V Cj ^\{j-l)r\ ci >0 V c. . (35) 

Then the Toeplitz matrix {M^;} = {^[(j —0^]} associated to the quadratic form (34) with arbitrary coefficients Cj 
has to be positive-definite. More explicitly, for even functions "tp, the Toeplitz matrix of order fc, 

/ V'(O) V’(^) '>p(2r) ... '>p[{k—l)r] 

ip{r) tp{0) . ip[{k — 2)r] 

\^/>[(fc—l)r] pj[{k—2)r] ... 'ip{r) ^/>(0) 

is positive-definite, with an eigenvalue spectrum bounded from below by zero^. 

As an instructive example of conditions resulting from the positive-definiteness of the matrices (36), let us consider 
the case of an even function pJ^r) and its corresponding 3x3 Toeplitz matrix, 

/ V'(O) -tpir) 'tp{2r)\ 

i){r) V’(O) ■0(0 • (37) 

y0(2r) 0(r) 0(0 )) 

Positive-definiteness implies positivity of the matrix determinant and of its minors along its diagonal. This gives the 
following set of inequalities, 

0(0) > 0(r), (38) 

0(0) > 0(2r) > ^|^-0(O), (39) 

where the last inequality comes from the determinant of (37), 

A = [0(0) — 0(2r)] [0^(0) — 20^(r) -|- 0(2r)0(O)] > 0 . (40) 

A practical method for positivity tests, coming from the straightforward generalization to higher order matrices, will 
be used in the following sections. 



The radial two-dimensional case 


Following an approach similar to the one-dimensional case, we want to relate the positivity (24) of the characteristic 
function to positive-definiteness properties of sets of matrices generalizing (but, actually, not of Toeplitz form) the 
matrices (36). 


® We are here using the polynomial method due to Riesz [8, 9]. On a more general footing, it comes from the application of a known 
theorem (see [9], Theorem 1.4): a necessary and sufficient condition that the trigonometric moment problem (31) have a generic (i.e. a 
solution whose spectrum is not reducible to a finite set of points) solution is that all Toeplitz quadratic forms 

k 

E fij-l Cj Cl > 0, k = 0,1,2^ ... , V complex vector c (33) 

3 , 1=0 

be positive. This means that the 2k x 2k matrices of moments {fj,j-i} are positive-definite, i.e. with smallest eigenvalue positive. 

^ The property (36) appears to be as a necessary consequence to the Bochner theorem [4] applied to the function '0(r) with a choice of 
points rj =j • r, j £ {1 k} Vk £ N Note that thanks to the r-dependence the condition (36) on the Toeplitz matrices ensures the 
positivity of the Fourier transform, as discussed in the previous section. 



Let us start with the coefficients of the Fourier series, 



1 



da dp 


F{a,P;r) , 


(41) 


recasting the characteristic function (24) with appropriate variables. For this sake, in analogy with the one-dimensional 
case, we consider real positive polynomials built from two complex variables z, z' and their complex conjugates, namely 


P{z,z\z,z') = \Q{z,z')f 


N 


2 


3,k 


E j ik—j'—ik' — ^ A 

Cj^k Z-'Z Z-’ Z Cji^k' > 0 . 

,k' 


(42) 


Choosing {z,z') = (e *“,e *^) and integrating (42) over F{a,p-,r), one finds. 




y 


n 2'K 

da dp F{a,p;r) 


y^ Cj,k ^ (r^y{j-j'Y + (fc-fc')^)) > 0 . 

j,k,j' ,k' 


(43) 

(44) 


Formula (44) implies the positive-definiteness of the tensorial form (44). A positive-definite matrix form can be 
obtained by noting that the polynomial Q{z, z') in (42) can be expanded [12] over the basis of monomials ordered by 
their degree. 






zz'2 , z'3 


(45) 


The positivity condition (44) reads as a positive, quadratic® form and thus as the definite-positiveness of an ordered 
hierarchy of matrices whose sizes depend on the chosen maximal orders of the corresponding polynomial (42). 

Let us illustrate this property by a low degree case, namely the positive definiteness of the 3x3 matrix, 

/ V’(O) i’ir) '4}{r\/2)\ 

I il){r) ■i/'(0) V’(^) I • (46) 

ip{r) iplo) ) 

Positive-definiteness implies positivity of the matrix determinant and of its minors along its diagonal, hence, 

-ipiO) > V’(r), (47) 

V'(O) > pJ{r^/2) > - V'(O) , (48) 

where the last inequality comes from the determinant of (46) 

A = [^(0) - [tjj{rV2)p;{0) - 2'4?{r) + V'^O)] > 0 . (49) 

Comparing with the similar one-dimensional case (40), it is worth noting that the minors’ inequalities from the matrix 
minors are the same as the one-dimensional ones (43), up to a rescaling of r. It is not the case for the last inequality 
of (48), coming from the determinant (49). Indeed, differences obviously occur because the new matrices, starting 
with (46) and beyond, are not any more of a Toeplitz type. Such differences will be the common rule at higher orders. 

The structure of the set of inequalities (47,48) can be elucidated by remarking that it stems from the application 
of the Bochner theorem to the set of positions, 


x) = {0,0} , {0,r} , {r,0}, (50) 

in a 2-dimensional square lattice. Indeed, using the 2-dimensional Bochner theorem [4], one finds the related necessary 
condition of positive-definiteness on the matrices (here a 3 x 3 matrix) of 2-vectors 

Wj,i} = . (51) 


® It is interesting to note that polynomials of two variables are not necessarily sums of squares but can always be expressed as a ratio of 
sum of squares [12]. Hence it is enough to ask for an arbitrary squared polynomial (42). 
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It is easy to recognize the identity of (51) with the matrix (46). Moreover typical minor’s inequalities correspond 
to the Bochner theorem for the pairs of points ({0, 0}, {1,0}) and ({0, 0}, {1,1}), respectively. This explains their 
relation with the one-dimensional properties. Putting together the three points {0, 0}, {1,0}, {1,1}, which are not 
aligned, gives rise to the matrix (46) which is not of Toeplitz form, characteristic of the one-dimensional problem. 

The same arguments easily extend to higher orders. For instance, at the next level, d = 2, one comes to a 6 x 6 
positive-definite matrix of 2-vectors {'ip{xj—xi)'\ corresponding to the basis (45) with six 2-vectors 

X- = {0,0} , {0,r} , {r,0} , {0,2r} , {r,r} , {2r,0} . (52) 

As in the one-dimensional case, the generalization to higher degrees is relatively straightforward and will lead to a 
subsequent application to Fourier-positivity in the radial two-dimensional case. 


V. APPLICATIONS TO FOURIER-POSITIVITY 

As theoretically motivated and developed in the preceding sections, Fourier-positivity of a real and even, positive 
function ^(r) can be tried and checked using in a finite set of rescalings of tp, namely S = {tp^kr)}, k = 1,..., K. In 
a first way that we call in short “Bochner method”, we make use of the positive-definiteness of r-dependent matrices 
whose components are given by S, as discussed in section IV, see (36) for the one-dimensional case and (51) for the 
radial two-dimensional case. 

The second way to test Fourier-positivity using the tools of section IV makes a direct use of the characteristic 
functions stemming from the Poisson resummation formulas, see respectively (10) for one-dimensional cases and (25) 
for the radial two-dimensional cases. The idea is to look for the appropriate range of the variable r, see resp. (20) 
and (29), for which the reconstruction of the characteristic functions from a finite set § satisfies the selection of the 
Fourier transform (p{s) with sufhcient accuracy to detect possible violations of positivity in some range of s. 

We shall test the capacity of the two different methods to thoroughly check Fourier-positivity or its violation. For 
this sake we introduce a large testing set of real, even and positive functions ipir), Fourier-positive or not, made of 
random combinations of a finite basis with well-known analytic Fourier transforms. To be more precise, we are using, 
in the one-dimensional case, the orthogonal basis of Hermite-Fourier functions, i.e. the quantum oscillator eigenstates 
which are eigenstates of the Fourier transform with eigenvalues 1 and -1. In the radial two-dimensional case we 
opt for an orthogonal basis of Laguerre polynomials multiplied by a simple exponential; these, under Fourier-Bessel 
transform, return combinations of rational functions fast decreasing at infinity. The random nature of the chosen 
combinations allow us to start with a large corpus of Fourier-positive and non Fourier-positive functions pj, allowing 
us to test our methods with a good accuracy and for quite different sets of Fourier partners. 


A. Fourier-Positivity in one dimension 

Starting with a large set of real even test functions, we examine the performance of the Fourier-positivity tests 
corresponding successively to the “Bochner method” and the “Poisson method”. 


A one-dimensional basis of Hermite-Fourier test functions 


Consider the Hermite-Fourier functions 

Up{r) = 7r"i e~i^^Hp{r). (53) 

Here, we set Hp to be a square normalized Hermite polynomial, with a positive coefficient for its highest power term. 
For the sake of clarity, we list the first polynomials as, Hq = I, Hi = ^/2 r, H 2 = {2r'^ — I)/-\/2, H 3 = (2r^ — 3r)/-\/3, 
and their recursion relation, 


^p+iHp-\-i — 2rapHp 2nap—iHp—i , 

where Op = y/2Pp\ . It is known that the Fourier transform of such states brings only a phase. 


(54) 


f dr Up{r) = Up{s ), 

\'27T J 


—00 


(55) 
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FIG. 1: Examples of “positive-positive” and “positive-negative” test functions. Left: both tppp (full line), with components, 
{.901, .276, .259, .006, .214}, and (p (dashed line) are positive. Right: ippn, with components, {.772, .304, .386, .171, .366}, remains 
> 0 but ip takes both signs. 


and thus such states give generalized self-dual functions with phase If one expands ip in the oscillator basis, 
'0(f) = '"P ^p(^)> a truncation at some degree N, then all odd order components C 2 p+i must vanish if p 

must be real, and the even rest splits, under Fourier transform, into an invariant part and a part with its sign reversed, 
namely 


[iV/4] lN/4] 

0(^) = X! “4p(f) + ^ C4p-2 U4p-2{r), 

p—0 p—1 

lN/4] lN/4] 

‘fis) = ^ C4p U4p{s) - C4p-2 U4p-2{s) , (56) 

p—0 p—1 

where the usual symbol [fV/4] means the integer part of N/A. This polynomial parametrization makes it trivial to 
generate fully positive ^/>’s, with both cases of partners p^s fully positive or p^s showing both signs. 

In our numerical illustration we consider the basis with N = 8, i.e. random combinations of the first five real eigen¬ 
states of the harmonic oscillator® whose coefficients, (cq, C 2 , C 4 , cq, cs}, are random real numbers with a normalization 
constraint, Cq -|- c| -b c® -b c§ -b c| = 1, and the condition that pj{r) > 0, Vr. We retained a randomly generated set of 
15456 positive functions ip{r) among which 4388 cases where p is also always positive and 11068 cases where p has a 
range of negative values. 


One-dimensional Bochner method 

For further discussions, let us denote 0pp(p), the “positive-positive” (resp. “positive-negative” ippn) test functions 
which lead to p positive (resp. not everywhere positive). Figure 1 shows a typical example in each category. 

Following the results of section IV, all Toeplitz matrices (36) built from 0pp(r) are positive definite, for all dimensions 
and for all values of the parameter r. None of our numerical tests based on positive-definiteness with our sample of 


® Our explicit parametrization is, 

tp{r) = 7 r“ 4 e“ 2 ’' [cq -t C 2 2 “ 2 ( 2 r^ — 1 )-|-C 4 6 “ 2 (4r'^ — 12 r^-|-3)/2 + 

C 6 5" ^ ( 8 r® - 60r4 -|- QOr^ - 15)/12 + cs 70" i (16r® - 224r® -|- SlOr^ - 840r2 + 105)/24 ] , (57) 

and with Fourier transform, 

(^(s) = 7r“ 4e“2® [co — C 2 2~2 (2s^ — 1) + C 4 6~2 (4s'^ — 12s^ + 3)/2 — 
ce5~i ( 8 s® - 605"^ + 90s^ - 15)/12 -I- cs 70" 5 (16s® - 224s® -I- 840s^ - 840s2 -I- 105)/24 ] . 


( 58 ) 







11 



FIG. 2: Lowest eigenvalue X as a function of the “grid” parameter r of the Toeplitz matrices, see (36). Upper curves, matrix 
dimension 5. Lower curves, matrix dimension 10. Left: tppp case, as described by left part of Fig.l. Right: 'ippn case shown in 
right part of Fig.l, non positivity of (p undetected by dimension 5, detected by dimension 10. 


4388 functions ippp{r) contradicted this fact. 

According to the same properties, the lack of positive definiteness of the Toeplitz matrix corresponding to the 
functions ippnif) will show off, sooner or later, for some dimension and some range of r. The sign of the Toeplitz 
determinant, however, may be misleading if an even number of negative eigenvalues occurs^*^. It is safer to track the 
lowest eigenvalue of a high enough dimensional matrix of the hierarchy (36). 

Figure 2 shows, for dimensions 5 and 10, how this lowest eigenvalue, A 5 ,Aio, respectively, behaves when r varies 
when choosing the example functions of Fig.l. Obviously, Aio < A 5 , since the 5-dimensional Toeplitz matrix is 
embedded in any higher-dimensional matrix of the hierarchy (36). The left part of Fig.2 considers again that case 
that was described by the left part of Fig.l. As should be, no negativity is observed, and a confirmation is expected 
for any matrix dimension. In turn, the right part of Fig.2 describes that case that was illustrated by the right part of 
Fig.l. No detection occurs if the Toeplitz matrix has only dimension 5, while dimension 10 provides a clear detection. 

Across our sample of 11068 cases, a proportion of ^ 69% are detected with Toeplitz matrices of dimension 5. 
With dimension 10, the success rate reaches ^ 86 %. The deeper the negative parts of </?, the higher the detection 
probability. However, the “Bochner method” may require for some “rebel” functions a very high dimensionality, and 
then becomes uneasy. The “Poisson method” will turn out to be easier and with almost full detection success. 


One-dimensional Poisson method 

Consider the parameter r as an integration grid parameter, Ar, and also the ratio, 0/Ar, as a pseudo momentum, 
s. Using the renormalized definition of the characteristic function (10) and the approximation (19) discussed in 
subsection 111 A, the range of tp in our samples allows a truncation into a finite sum, namely, 

K 

F(s,Ar) ~ ArjV^ ip{nAr) exp[f(nAr)s] (59) 

n=-K 

with K ~ R/{Ar). Here a range i? = 10 is enough to perform a correct coverage of the characteristic function 
F{s,Ar). Note that we thus keep KAr ~ R constant and look for an integration grid parameter Ar not too small, 
bounded from below, as discussed in subsection HI A. 

Contour lines of the values of F, in terms of Ar and s, are shown in Figure 3. Its left and right parts correspond to 
the “doubly positive” and “partly negative” cases already used for the previous Figures. As expected, no contour for 
a negative value of F is found for the “doubly positive” case, while “negative contours” occur for the case when if is 
partly negative. Note, for instance, the “negative” contour line where F = —0.028. A few big dots have been plotted 
to reinforce the visual identification of such “negative” contours. 


A way to evade this difSculty would be to increase the dimension by one, but this costing computation time. 
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FIG. 3: Contour plots of the characteristic function F in terms of the discretization parameter Ar and the pseudo-momentum 
s = Of Ar. Left: Double positivity case, already described by left part of Fig.l. No contour line is found for negative values of 
F. Right: Case already shown in right part of Fig.l, non positivity of ip detected by contours for negative values of F. 


Such a result, namely a good approximation to a brute force Fourier transform, is expected when Ar is small enough 
to ensure a good convergence, F’(s, Ar) —>■ ip{s), but larger values of Ar maintain the criterion: negative values of F 
occur only for non positive (/?’s. 

A comment is in thus order: It can be remarked from Fig.3 that the contour curves remain parallel to the abscissa 
axis for a rather large range of Ar. This, joined to the previous remark, can be explained by the fact that the Fourier 
transform ip{s) can be quite well reconstructed from the approximation (59) independently from the value of Ar or 
equivalently of the summation number K. This is exemplified in Figure 4 for the non Fourier-positive function of Fig.l 
(right). This is why negative values are reproduced and the test for non Fourier-positivity is successful. 

We verified, for Ar running between .1 and 1 and 0 < s < 8 , that our 4388 “doubly positive” test cases did not 
generate any negative value of F. In turn we verified, for Ar running between .1 and 1 and 0 < s < 8 again, that most 
partly negative (p's are detected by negative contours. Typically, for 11068 test cases, only 19 of them fail generating 
negative values of F and thus escape detection. An inspection of such “rebel” cases gives a clear explanation for the 
failure: the negative values of such p’s are tiny. 

It can be concluded that F provides a very efficient test for the selection of ip without a detailed calculation of p. 


B. Fourier-Bessel-positivity 

Basis of radial funetions conneeted by Bessel transform 

Using a similar method as in one dimension, we set a basis, for the functions 'll), built from the exponential, e~^, 
multiplied by random linear combinations of nine Laguerre polynomials^^. These are normalized, x'il!{x)‘^dx = 1, 
from the condition, cf = 1, for the random, real number, mixture coefficients, Cj, i = 0,..., 8. The same coefficients 
are also selected so that ^(x) be positive, Vx > 0. The partners in radial momentum space, p{p) = xJo(px)'ip(x)dx, 


The explicit form of the basis is 

e“i x{l, {-2 + x)lV2, (io-Qx + :F)/{2y/Z), (-24-h 36a: - 12x2-|-x3)/12, (i20-240x-P120x^-20x^-Px^)/(24V5), (-720-|-1800x- 
1200x2-1-300x3-30x"‘-|-x®)/(120'/6), (5040-15120:^+12600x2-4200x3-|-630x"‘-42x®-|-x®)/(720U7), (-40320+141120x-141120x2 + 
58800x3 _ 11700x4 + 1176x® - 56x® + x'^)/(10080V2), (362880 - 1451520x + 1693440x2 - 846720x3 + 211680x4 - 28224x5 + 2016x6 - 
72x7 + x«)/120960)}. 
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FIG. 4: Example of a ‘reconstruction” of a non-Fourier-positive test function. The function (fpn, depicted in the right part of 
Fig.l by a dashed line, is approximated respectively with K = 12 (small dashes), IF = 14 (large dashes) and IF = 20 (full line, 
indistinguishable from the true ippn-)- 


are obtained, with the same coefficients Ci, from the corresponding basis of functions^^ , also normalized. It is clear 
that i.p{p) reads as a polynomial divided by a common denominator, (1 + This makes it easy to sort out 

positive Lp's from those which take both negative and positive values. We show in Figure 5 a double positivity case 
(left part) and a case with Lp partly negative (right part). 

Our test basis for double positivity contains 185 cases and that for situations where only if remains always positive 
contains 9894 cases. 








FIG. 5: Radial partners under Bessel transform, if (full lines) and ip (dashes). Left : Fourier-positive case if = e 2 (3.6096 — 
6.7462a: -f 4.8826®^ - 1.6141x® -f 0.28086a;‘‘ - 0.027009a;® -b 0.00145a;® - 0.000042169x^ -b 5.63539 10"’^a;®). Right : non Fourier¬ 
positive case ip = e-'i (2.49362 - 6.84573a:-b6.76697x^ -3.04127x® -b0.723816a:‘‘ -0.0959944x® -b0.00705616x® -0.000265057x’' -b 
3.93896 10~®a:®). For graphical reasons, the Figure actually shows ;/;(p/7)/7. 


The basis in Fourier transformed space reads : 

{4/(1 + 4p2)3/2^ (_4^(_1 + 8p2))/(l -I- 4p2)5/2^ (4v/3(l - 24p2 + 48p'*))/(l -b 4p2)V2^ (-8(-l -b 48p2 - 288p'‘ + 256p®))/(l -b 
4p2)9/2^ (4^(1 _ 80p2 -bQeOp-* - 2560p® -b I280p®))/(1 -b 4p2)ii/2^ + i20p2 - 2400p‘‘ -b 12800p® - 19200p® -b 6144pi0))/(l -b 

4 p2jl3/2_ (4v^(l - 168p2 + 5040p‘‘ - 44800p® -b 134400p® - 129024pl° + 28672pl2))/(l + 4p2)l5/2, (_8y2(-l -b 224p2 - 9408p‘‘ -b 
125440p® - 627200p® -b 1204224pi° - 802816pi2 -b 131072pi‘‘))/(i + 4p2)i7/2^ ( 42(1 - 288p2 -b 16128p‘‘ - 301056p® -b 2257920p® - 
7225344pio + 9633792pi2 - 4718592pi4 -b 589824pi®))/(l -b 4p2)i9/2}, 















14 




FIG. 6: Evolution of the lowest eigenvalue of a Bochner matrix as a function of the scaling parameter /3 for the random points 
which parametrize the matrix. From top to bottom curve, dimension 20,40,60,80. Left: the eigenvalue soon becomes negative 
and detects a partly negative ip. Right: failure case, where a much larger dimension would be needed for the detection. 


Bochner method for radial functions 

It must be kept in mind here that we are in a two-dimensional situation, namely that, given a point with coordinates 
{y,z}, the argument a; of '0 is, x = \/y^ + z?. Given a list of points {yi,Zi},i = 1,..., FT, the matrix elements of 
the associated, iLth-order Bochner matrix read, 0 (y^{yi — yjY + (zi — Zj)^^ . Our numerical tests used a set of points 

fi the coordinates of which, {yi,Zi}, are random real numbers in a range, {—20,20}. A scale parameter, /3, is then 
introduced to adjust such points to any range {—20/3,20/3}. Typically, we considered the first 20 points, {/3yi,/3zi}, 
with for instance (3 = .5, but we also used /3 = .1, 0 = .4, /3 = .8, f3 = 1. When 20 points gave too few detections, we 
used as many as 80 or even 100 points. 

As expected, the Bochner matrices are found positive definite when we investigate the 185 “doubly positive” cases. 
In turn, our set of 9884 test functions for partly negative Lp's returned a detection rate of ^ 20% when Bochner 
matrices of dimension 20 were used. The rate reached ^ 40% for matrices of dimension 80 and hardly increased if 100 
points were used. This modest detection rate with reasonable size matrices contrasts with the good result obtained 
in the one-dimensional, Gaussian test function case. Matrices of a much higher order are now needed, or a set of 
more efficient points r} must be defined. But we tried several sets of points, different from random ones, and failed to 
design a “maximum efficiency set of points”. 

Fig. 6 shows, for dimensions 20,40, 60,80 respectively, the evolution of the lowest eigenvalue A of a Bochner matrix 
as a a function of the scale parameter /3. Naturally, the matrix with dimension 20 being a submatrix of that with 
dimension 40, the latter generates a lower bound. The same reasoning applies when the dimension increases, hence 
the ordering of the four shown curves. The left part of the Figure describes a case where negative values of A are fast 
obtained. The left part describes a failure case (even for dimension 100). 


Poisson method for radial functions 

With the obvious symmetries at our disposal, the Poisson function in this situation can be rewritten as, 

K K _ 

F{a,"f,Ar) = (Ar)^/(27r) EE emen'fp (^Ar\/m? + cos(toq:) cos{n'y) , (60) 

m—0 k—0 

where, = 2 — Smo, and the same for e„ , account for edge effects. The range i? of 0, typically i? ~ 40 for our 
test functions, provides a natural cut-off, K ~ R/Ar, for the {m,n} summations. The normalization coefficient, 
(Ar)^/(27r), has been introduced to make F similar to a Fourier integral at the limit, Ar —^ 0. 

A similar comment to the one-dimensional case is in order for the radial case. An approximate reconstruction of 
the Fourier transform function appears to be allowed following formula (19) 

Figure 7 shows that, for a finite value of Ar, hence for finite summations governed by the resulting K, negative 
values of F can be found, even for a “rebel case” like that shown in the right part of Fig. 6. A search for negative 
values of F under moderate values of Ar returns a detection rate of ~ 90% through our set of 9884 functions. This is 
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FIG. 7: Contour plot of F{a,y) when Ar = .5. Left: same case as the “success” case shown in the left part of Fig.6; contours 
for F = —.3,3 are shown. Right: same case as the “failure” case shown in the right part of Fig.6; but now, a contour with a 
negative value, F = —.02, provides a detection. 





FIG. 8: Example of a “reconstruction” of a radial non-Fourier-positive test function. The function ippn depicted in the right 
part of Fig.5 by a dashed line is approximated respectively with K = 40 (large dashes, detecting one negative domain) and 
with K = 80 (full line, detecting both negative domains). 


significantly better than the result discussed in the previous subsection, where the practical tests based on “Bochner 
method” appears to be only very slowly evolving with already high matrix order. 

Fig.7 calls for a comment. The plots appear to approximately satisfy a radial symmetry in the two-dimensional 
(a, 7 ) plane, while the resummation formula (60) is not a priori symmetric. The reason is that it gives an approximate 
reconstruction of the radial function ’p{s) as predicted by our general argument of section II for the radial 2-dimensional 
case. This reconstruction property with a finite number of terms in (60) is exemplified in Fig .8 with the non Fourier¬ 
positive function depicted in Fig.5. The negative domain is detected already with K = AQ and reasonably fully 
reproduced with K = 80. 


VI. SUMMARY AND OUTLOOK 

In any dimension d there exists a set of real functions, partners under Fourier transform, that are both positive. 
We have called them Fourier-positive functions. If only for purely mathematical curiosity, the mapping between these 
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two convex, partner sets is of some interest. But this interest is reinforced by the fact that, in theoretical physics, it 
happens that a “density” in one space has a Fourier partner which is also a “density”. (By definition, densities are 
positive observables.) 

There is also an interesting “subproblem”, that of the mapping between subsets, for instance polynomials multiplied 
by Gaussians in both spaces, or polynomials multiplied by simple exponentials in one space with rational function 
partners. Such subsets are nested according to the order of the considered polynomials. This hierarchy allows a useful 
set of successive approximations. 

In fact a general, both necessary and sufficient, mathematical criterion for ensuring Fourier-positivity seems not to 
be yet known. Indeed, examples taken from physics show that Fourier-positivity is a nontrivial constraint on density 
models [14], where small modifications may play a role. For a simple illustration, compare a square well density, 
whose Fourier transform is the Bessel function, with a Gaussian density, which, being invariant by Fourier transform, 
is obviously Fourier positive. 

Modern computers allow “Fast Fourier Transforms” which give an easy answer to positivity properties in both 
partner spaces, but, obviously, a pure numerical approach is not completely satisfactory. The present work gives 
several mathematical, analytical arguments to complement our previous work [5], where it was shown that the topology 
of such interesting “positive partner subsets” was highly non trivial and moreover, where there appeared an intuition 
that extremal elements in such convex sets are reminiscent of Dirac combs. 

The proofs displayed in this work do take advantage of this intuition, but indirectly. In substance, we use two kinds 
of criteria to test whether a positive “object” ip has a positive “image” if, criteria that use values of ip only, i) A 
positivy criterion for the characteristic Poisson function associated with ip through a “Dirac comb” distribution turns 
out to be quite efficient for both d = 1 and radial, d = 2 cases, ii) A criterion taken from Boechner’s theorem, namely 
the positivity of Toeplitz matrices and similar matrices, turns out to be efficient for d = 1 cases, but disappointing for 
radial, d=2 cases. In both cases, we took great care to validate our analytical considerations by means of controled, 
numerical tests, including statistical evaluations. 

A third set of criteria [6] was in fact our initial approach. It consisted in relating positivity of a function and 
convexity of analytical continuations of related functions, to define bounds via the Jensen’s theorem, but we failed 
in making this approach a convincing one, if only because analytical continuation most often has to face severe 
singularities. It is not displayed in the present work. We keep it on a back-burner. 

On a deeper but difficult level, the question of a general criterion of Fourier positivity and a classification of those 
functions still remains widely open. We hope that the consideration of Dirac combs and Poisson resummation proposed 
in our paper in this context may help to make some new steps in that problem. 

For a more immediate outlook, our priority will be to take advantage of the hierarchy of mapped subspaces above- 
mentioned and use it for actual physical problems. Reliable error bars are essential in data analysis and we want to 
use our criteria for both estimates of ip and ip. 
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